library(ggplot2)
library(haven)
library(plyr)
library(dplyr)
library(ggpubr)
library(gridExtra)
library(grid)
library(stargazer)
library(lme4)
library(optimx)
library(lmerTest)
library(nlme)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(broom.mixed)
library(dotwhisker)
library("psych")
library(merTools)
library(ggeffects)


#Experimental analsysis, wide-format data
load("replication_data_UA.rda")

#Descriptive results
data_msd <- WIDECJ %>%                           # Get mean & standard deviation by group
  group_by(F3, F2, F1) %>%
  summarise_at(vars(RATING),
               list(mean = mean,
                    sd = sd,
                    se = ~sd(.)/sqrt(n())), na.rm = TRUE) %>% 
  as.data.frame()

data_msd  

data_msd$F2 <- factor(data_msd$F2, levels = c("US supply", "US withdraw"))
data_msd$F1 <- factor(data_msd$F1, levels = c("continuing", "intensifying", "extending non-EU","extending EU"))
data_msd$F3 <- factor(data_msd$F3, levels = c("maintain", "increase independently","increase cooperation","integrated army"))


ggplot(data_msd,                               
       aes(x = F1, shape = F2,
           group = F3,
           y = mean)) + 
  geom_errorbar(aes(ymin=mean-2*se, ymax=mean+2*se), width = 0.1)+
  facet_wrap(~F3)+
  geom_point(size = 5)+
  theme_minimal()+  
  theme(axis.text.x=element_text(angle=30,hjust=1), text = element_text(size = 17))


WIDECJ$F1 <- factor(WIDECJ$F1, levels = c("continuing", "intensifying", "extending non-EU","extending EU"))
WIDECJ$F2 <- factor(WIDECJ$F2, levels = c("US supply", "US withdraw"))
WIDECJ$F3 <- factor(WIDECJ$F3, levels = c("maintain", "increase independently","increase cooperation","integrated army"))

mod0 <- lmer(RATING ~ F1 + F2 + F3 + (1|ID) + as.factor(Country), data = WIDECJ, REML=TRUE,
             control = lmerControl(
               optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(mod0)
anova(mod0)

mod1 <- lmer(RATING ~ F1*F2*F3 + (1|ID) + as.factor(Country), data = WIDECJ, REML=TRUE,
             control = lmerControl(
               optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(mod1)
anova(mod1)

#stargazer(mod1, out = "3waymodel.html",type = "html")
class(mod0) <- "lmerMod"
stargazer(mod0, out = "simplemodel.tex",type = "latex")

m1_df <- tidy(mod1)

four_brackets <- list(
  c("Main Effects", "F1intensifying", "F3integrated army"),
  c("Country FEs", "as.factor(Country)France", "as.factor(Country)Portugal"),
  c("2-way*", "F1intensifying:F2US withdraw", "F2US withdraw:F3integrated army"),
  c("3-way*", "F1intensifying:F2US withdraw:F3increase independently", "F1extending EU:F2US withdraw:F3integrated army")
)

two_brackets <- list(
  c("Main Effects", "F1intensifying", "F3integrated army"),
  c("Country FEs", "as.factor(Country)France", "as.factor(Country)Portugal")
)


p1 <- {dwplot(list(first=mod1), effects="fixed")+
    geom_vline(xintercept=0, lty=2)+
    scale_y_discrete(labels=c("as.factor(Country)Germany" = "Germany", 
                              "as.factor(Country)Hungary" = "Hungary",
                              "as.factor(Country)Italy" ="Italy",
                              "as.factor(Country)France" ="France",
                              "as.factor(Country)Poland" ="Poland",
                              "as.factor(Country)Portugal" ="Portugal"))+
    scale_colour_grey()+theme_minimal()+theme(legend.position = "none")} %>%
  add_brackets(four_brackets, fontSize = 0.65)

p2<-{dwplot(list(first=mod0), effects="fixed")+
    geom_vline(xintercept=0, lty=2)+
    scale_y_discrete(labels=c("as.factor(Country)Germany" = "Germany", 
                              "as.factor(Country)Hungary" = "Hungary",
                              "as.factor(Country)Italy" ="Italy",
                              "as.factor(Country)France" ="France",
                              "as.factor(Country)Poland" ="Poland",
                              "as.factor(Country)Portugal" ="Portugal"))+
    scale_colour_grey()+theme_minimal()+theme(legend.position = "none")} %>%
  add_brackets(two_brackets, fontSize = 0.65)

p1

p2

# Interaction F1 and F3:

smesF1 <- lapply(split(WIDECJ,WIDECJ$F1),
                 lmer,formula=RATING ~ F3 + (1|ID)+ as.factor(Country))

smesF1

p1 <- dwplot(list(first=smesF1[[1]]), effects="fixed")+
  geom_vline(xintercept=0, lty=2)+
  ggtitle("continuing")+ 
  ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
  xlim(-0.5,1.5)+
  scale_colour_grey()+theme_minimal()+theme(legend.position = "none")
p1
p2 <- dwplot(list(first=smesF1[[2]]), effects="fixed")+
  geom_vline(xintercept=0, lty=2)+
  ggtitle("intensifying")+ 
  ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
  xlim(-0.5,1.5)+
  scale_colour_grey()+theme_minimal()+theme(legend.position = "none")
p2
p3 <- dwplot(list(first=smesF1[[3]]), effects="fixed")+
  geom_vline(xintercept=0, lty=2)+
  ggtitle("extending non-EU")+ 
  ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
  xlim(-0.5,1.5)+
  scale_colour_grey()+theme_minimal()+theme(legend.position = "none")
p3
p4 <- dwplot(list(first=smesF1[[4]]), effects="fixed")+
  geom_vline(xintercept=0, lty=2)+
  ggtitle("extending EU")+ 
  ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
  xlim(-0.5,1.5)+
  scale_colour_grey()+theme_minimal()+theme(legend.position = "none")
p4
grid.arrange(p1,p2,p3,p4)


class(smesF1$continuing) <- "lmerMod"
class(smesF1$intensifying) <- "lmerMod"
class(smesF1$`extending non-EU`) <- "lmerMod"
class(smesF1$`extending EU`) <- "lmerMod"

stargazer(smesF1$continuing, out = "F1F3continuing",type = "latex")
stargazer(smesF1$continuing, out = "F1F3continuing",type = "latex")
stargazer(smesF1$continuing, out = "F1F3continuing",type = "latex")
stargazer(smesF1$continuing, out = "F1F3continuing",type = "latex")


dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")
dat4 <- ggemmeans(smesF1[[4]], terms = c("F3"), type="fixed")



df <- bind_rows(
  dat1 %>% mutate(model = "continuing"),
  dat2 %>% mutate(model="intensifying"),
  dat3 %>% mutate(model="extending non-EU"),
  dat4 %>% mutate(model="extending EU")
)

pp1 <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model), position=position_dodge(width=0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model), position=position_dodge(width=0.2))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("continuing"="grey0","intensifying"="grey30","extending non-EU"="grey50","extending EU"="grey80"), name = "F1 - Russia actions")+
  coord_flip()

pp1

df <- df[,c(1:5,7)]
stargazer(df, summary = FALSE, out = "F1F3predicted",type = "latex")

# Interaction F2 and F3:

smesF1 <- lapply(split(WIDECJ,WIDECJ$F2),
                 lmer,formula=RATING ~ F3 + (1|ID)+ as.factor(Country))

smesF1

p1 <- dwplot(list(first=smesF1[[1]]), effects="fixed")+
  geom_vline(xintercept=0, lty=2)+
  ggtitle("US supply")+ 
  ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
  xlim(-0.5,1.5)+
  scale_colour_grey()+theme_minimal()+theme(legend.position = "none")
p1
p2 <- dwplot(list(first=smesF1[[2]]), effects="fixed")+
  geom_vline(xintercept=0, lty=2)+
  ggtitle("US withdraw")+ 
  ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
  xlim(-0.5,1.5)+
  scale_colour_grey()+theme_minimal()+theme(legend.position = "none")
p2

grid.arrange(p1,p2)


class(smesF1$`US supply`) <- "lmerMod"
class(smesF1$`US withdraw`) <- "lmerMod"

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

pp1 <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2 - US reaction")+
  coord_flip()

pp1

df <- df[,c(1:5,7)]
stargazer(df, summary = FALSE, out = "F2F3predicted",type = "latex")


# Interaction F2 and F3 by country:
for (i in 1:length(names(table(WIDECJ$Country)))){
  COUNTRYX <- WIDECJ[WIDECJ$Country %in% names(table(WIDECJ$Country))[i],]
  
  smesF1 <- lapply(split(COUNTRYX,COUNTRYX$F2),
                   lmer,formula=RATING ~ F3 + (1|ID))
  
  smesF1
  
  
  p1 <- dwplot(list(first=smesF1[[1]]), effects="fixed")+
    geom_vline(xintercept=0, lty=2)+
    ggtitle("US supply")+ 
    ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
    xlim(-1,2.5)+
    scale_colour_grey()
  p1
  p2 <- dwplot(list(first=smesF1[[2]]), effects="fixed")+
    geom_vline(xintercept=0, lty=2)+
    ggtitle("US withdraw")+ 
    ylim(breaks=c("F3increase independently","F3increase cooperation","F3integrated army"))+
    xlim(-1,2.5)+
    scale_colour_grey()
  p2
  grid.arrange(p1,p2, top = names(table(WIDECJ$Country))[i])
}

# Interaction ideology and F3:

smesF1 <- lapply(split(WIDECJ,WIDECJ$ideology_rec2),
                 lmer,formula=RATING ~ F3 + (1|ID)+ as.factor(Country))

smesF1

class(smesF1$left) <- "lmerMod"
class(smesF1$center) <- "lmerMod"
class(smesF1$right) <- "lmerMod"

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "left"),
  dat2 %>% mutate(model="center"),
  dat3 %>% mutate(model="right")
)

pp1 <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model), position=position_dodge(width=0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model), position=position_dodge(width=0.2))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("left"="grey0","center"="grey40","right"="grey70"), name = "Ideology")

pp1

df <- df[,c(1:5,7)]
stargazer(df, summary = FALSE, out = "F3ideology",type = "latex")

# Interaction veto and F3: ####

WIDECJ$veto <- WIDECJ$B13_1

WIDECJ$veto[WIDECJ$veto %in% c(99)] <- NA
WIDECJ$veto[WIDECJ$veto<=3] <- 1
WIDECJ$veto[WIDECJ$veto>3 & WIDECJ$veto<7] <- 2
WIDECJ$veto[WIDECJ$veto>=7] <- 3
WIDECJ$veto <- revalue(as.factor(WIDECJ$veto), c("1"="antiveto","2"="neutral", "3"="proveto"))

WIDECJ$veto 

smesF1 <- lapply(split(WIDECJ,WIDECJ$veto),
                 lmer,formula=RATING ~ F3 + (1|ID)+ as.factor(Country))

smesF1

class(smesF1$antiveto) <- "lmerMod"
class(smesF1$neutral) <- "lmerMod"
class(smesF1$proveto) <- "lmerMod"

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "antiveto"),
  dat2 %>% mutate(model="neutral"),
  dat3 %>% mutate(model="proveto")
)

pp1 <- ggplot(na.omit(df), aes(x, predicted, group=model)) +
  geom_point(aes(colour=model), position=position_dodge(width=0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model), position=position_dodge(width=0.2))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("antiveto"="grey0","neutral"="grey40","proveto"="grey70"), name = "Veto")

pp1

df <- df[,c(1:5,7)]
stargazer(df, summary = FALSE, out = "F3ideology",type = "latex")


# Interaction identity and F3:


smesF1 <- lapply(split(WIDECJ,WIDECJ$EU_identity_rec2),
                 lmer,formula=RATING ~ F3 + (1|ID)+ as.factor(Country))

smesF1

class(smesF1$`EU only`) <- "lmerMod"
class(smesF1$`EU & Nat.`) <- "lmerMod"
class(smesF1$`Nat. & EU`) <- "lmerMod"
class(smesF1$`Nat. only`) <- "lmerMod"

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")
dat4 <- ggemmeans(smesF1[[4]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "EU only"),
  dat2 %>% mutate(model="EU & Nat."),
  dat3 %>% mutate(model="Nat. & EU"),
  dat4 %>% mutate(model="Nat. only")
)

pp2 <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model), position=position_dodge(width=0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model), position=position_dodge(width=0.2))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("EU only"="grey0","EU & Nat."="grey40","Nat. & EU"="grey65","Nat. only"="grey75"), name = "Identity")

pp2

df <- df[,c(1:5,7)]
stargazer(df, summary = FALSE, out = "F3identity",type = "latex")



grid.arrange(pp1,pp2, ncol=2)

# Interaction threat and F1:


smesF1 <- lapply(split(WIDECJ,WIDECJ$threat_all_rec),
                 lmer,formula=RATING ~ F1 + (1|ID)+ as.factor(Country))

smesF1

class(smesF1$`no/small`) <- "lmerMod"
class(smesF1$`moderate`) <- "lmerMod"
class(smesF1$`large`) <- "lmerMod"

dat1 <- ggemmeans(smesF1[[1]], terms = c("F1"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F1"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F1"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "no/small"),
  dat2 %>% mutate(model="moderate"),
  dat3 %>% mutate(model="large")
)

pp2 <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model), position=position_dodge(width=0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model), position=position_dodge(width=0.2))+
  theme_minimal()+xlab("F1 - Russia actions")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("no/small"="grey0","moderate"="grey40","large"="grey65"), name = "Threat")

pp2


smesF1 <- lapply(split(WIDECJ,WIDECJ$threat_socio),
                 lmer,formula=RATING ~ F1 + (1|ID)+ as.factor(Country))

smesF1

class(smesF1$`no/small`) <- "lmerMod"
class(smesF1$`moderate`) <- "lmerMod"
class(smesF1$`large`) <- "lmerMod"

dat1 <- ggemmeans(smesF1[[1]], terms = c("F1"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F1"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F1"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "no/small"),
  dat2 %>% mutate(model="moderate"),
  dat3 %>% mutate(model="large")
)

pp2 <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model), position=position_dodge(width=0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model), position=position_dodge(width=0.2))+
  theme_minimal()+xlab("F1 - Russia actions")+ylab("Predicted Rating (0-10 scale)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("no/small"="grey0","moderate"="grey40","large"="grey65"), name = "Threat")

pp2




# Interaction F2XF3 by NATO trust:

WIDECJ$trustNATO <- WIDECJ$A4_3
WIDECJ$trustNATO <- as.numeric(WIDECJ$trustNATO)
WIDECJ$trustNATO[WIDECJ$trustNATO %in% c(98,99)] <- NA

WIDECJ$trustNATO_rec <- WIDECJ$trustNATO
WIDECJ$trustNATO_rec[WIDECJ$trustNATO_rec<=3] <- 1
WIDECJ$trustNATO_rec[WIDECJ$trustNATO_rec>3 & WIDECJ$trustNATO_rec<7] <- 2
WIDECJ$trustNATO_rec[WIDECJ$trustNATO_rec>=7] <- 3

WIDECJ$trustNATO_rec <- as.factor(WIDECJ$trustNATO_rec)
WIDECJ$trustNATO_rec <- revalue(WIDECJ$trustNATO_rec, c("1"="low", "2"="moderate", "3"="high"))


NATOX <- WIDECJ[WIDECJ$trustNATO_rec %in% names(table(WIDECJ$trustNATO_rec))[1],]

smesF1 <- lapply(split(NATOX,NATOX$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

low <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("NATO Trust",paste(names(table(WIDECJ$trustNATO_rec))[1]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()

low



NATOX <- WIDECJ[WIDECJ$trustNATO_rec %in% names(table(WIDECJ$trustNATO_rec))[2],]

smesF1 <- lapply(split(NATOX,NATOX$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

moderate <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("NATO Trust",paste(names(table(WIDECJ$trustNATO_rec))[2]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()

moderate


NATOX <- WIDECJ[WIDECJ$trustNATO_rec %in% names(table(WIDECJ$trustNATO_rec))[3],]

smesF1 <- lapply(split(NATOX,NATOX$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

high <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("NATO Trust",paste(names(table(WIDECJ$trustNATO_rec))[3]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()

high


grid.arrange(low, moderate, high, ncol=1)


# Interaction F2XF3 by pacifism:

table(WIDECJ$weapons)

i=1
ESCAL <- WIDECJ[WIDECJ$weapons %in% names(table(WIDECJ$weapons))[i],]

smesF1 <- lapply(split(ESCAL,ESCAL$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

good <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("Arms to Ukraine do mostly",paste(names(table(WIDECJ$weapons))[i]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()

good

i=2
ESCAL <- WIDECJ[WIDECJ$weapons %in% names(table(WIDECJ$weapons))[i],]

smesF1 <- lapply(split(ESCAL,ESCAL$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

neutral <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("Arms to Ukraine do mostly",paste(names(table(WIDECJ$weapons))[i]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()

neutral 

i=3
ESCAL <- WIDECJ[WIDECJ$weapons %in% names(table(WIDECJ$weapons))[i],]

smesF1 <- lapply(split(ESCAL,ESCAL$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

harm <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("Arms to Ukraine do mostly",paste(names(table(WIDECJ$weapons))[i]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()



grid.arrange(good, neutral, harm, ncol=1)


# Interaction F2XF3 by pacifism escalate:

table(WIDECJ$escalate)

ESCAL <- WIDECJ[WIDECJ$escalate %in% names(table(WIDECJ$escalate))[1],]

smesF1 <- lapply(split(ESCAL,ESCAL$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

no <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("join US efforts to weaken Russia, even at risk of escalating",paste(names(table(WIDECJ$escalate))[1]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


ESCAL <- WIDECJ[WIDECJ$escalate %in% names(table(WIDECJ$escalate))[3],]

smesF1 <- lapply(split(ESCAL,ESCAL$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

neutral <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("join US efforts to weaken Russia, even at risk of escalating",paste(names(table(WIDECJ$escalate))[3]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


ESCAL <- WIDECJ[WIDECJ$escalate %in% names(table(WIDECJ$escalate))[2],]

smesF1 <- lapply(split(ESCAL,ESCAL$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

yes <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("join US efforts to weaken Russia, even at risk of escalating",paste(names(table(WIDECJ$escalate))[2]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


grid.arrange(no, neutral, yes, ncol=1)


# Interaction F1XF3 by pacifism:

table(WIDECJ$weapons)

ESCAL <- WIDECJ[WIDECJ$weapons %in% names(table(WIDECJ$weapons))[1],]

smesF1 <- lapply(split(ESCAL,ESCAL$F1),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")
dat4 <- ggemmeans(smesF1[[4]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "continuing"),
  dat2 %>% mutate(model="intensifying"),
  dat3 %>% mutate(model="extending non-EU"),
  dat4 %>% mutate(model="extending EU")
)

good <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("continuing"="grey0","intensifying"="grey30","extending non-EU"="grey50","extending EU"="grey80"), name = "F1")+
  ggtitle("Arms to Ukraine do mostly",paste(names(table(WIDECJ$weapons))[1]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


ESCAL <- WIDECJ[WIDECJ$weapons %in% names(table(WIDECJ$weapons))[2],]

smesF1 <- lapply(split(ESCAL,ESCAL$F1),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")
dat4 <- ggemmeans(smesF1[[4]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "continuing"),
  dat2 %>% mutate(model="intensifying"),
  dat3 %>% mutate(model="extending non-EU"),
  dat4 %>% mutate(model="extending EU")
)

neutral <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("continuing"="grey0","intensifying"="grey30","extending non-EU"="grey50","extending EU"="grey80"), name = "F1")+
  ggtitle("Arms to Ukraine do mostly",paste(names(table(WIDECJ$weapons))[2]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()



ESCAL <- WIDECJ[WIDECJ$weapons %in% names(table(WIDECJ$weapons))[3],]

smesF1 <- lapply(split(ESCAL,ESCAL$F1),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")
dat3 <- ggemmeans(smesF1[[3]], terms = c("F3"), type="fixed")
dat4 <- ggemmeans(smesF1[[4]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "continuing"),
  dat2 %>% mutate(model="intensifying"),
  dat3 %>% mutate(model="extending non-EU"),
  dat4 %>% mutate(model="extending EU")
)

harm <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("continuing"="grey0","intensifying"="grey30","extending non-EU"="grey50","extending EU"="grey80"), name = "F1")+
  ggtitle("Arms to Ukraine do mostly",paste(names(table(WIDECJ$weapons))[3]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()

grid.arrange(good, neutral, harm, ncol=1)

# Interaction F2XF3 by threat:

table(WIDECJ$threat_all_rec)

i=3
THREAT <- WIDECJ[WIDECJ$threat_all_rec %in% names(table(WIDECJ$threat_all_rec))[i],]

smesF1 <- lapply(split(THREAT,THREAT$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

large <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("Threat",paste(names(table(WIDECJ$threat_all_rec))[i]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


i=1
THREAT <- WIDECJ[WIDECJ$threat_all_rec %in% names(table(WIDECJ$threat_all_rec))[i],]

smesF1 <- lapply(split(THREAT,THREAT$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

no <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("Threat",paste(names(table(WIDECJ$threat_all_rec))[i]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


i=2
THREAT <- WIDECJ[WIDECJ$threat_all_rec %in% names(table(WIDECJ$threat_all_rec))[i],]

smesF1 <- lapply(split(THREAT,THREAT$F2),
                 lmer,formula=RATING ~ F3 + (1|ID))
smesF1

dat1 <- ggemmeans(smesF1[[1]], terms = c("F3"), type="fixed")
dat2 <- ggemmeans(smesF1[[2]], terms = c("F3"), type="fixed")

df <- bind_rows(
  dat1 %>% mutate(model = "US supply"),
  dat2 %>% mutate(model="US withdraw")
)

moderate <- ggplot(df, aes(x, predicted, group=model)) +
  geom_point(aes(colour=model)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1, color = model))+
  theme_minimal()+xlab("F3 - capacity")+ylab("Pred.Rating (0-10)")+
  theme(legend.position = "right")+
  scale_color_manual(values=c("US supply"="grey0","US withdraw"="grey50"), name = "F2")+
  ggtitle("Threat",paste(names(table(WIDECJ$threat_all_rec))[i]))+ 
  theme(axis.text.x = element_text(angle = 45))+ylim(3.5,8.5)+
  coord_flip()


grid.arrange(no, moderate, large, ncol=1)

